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ABSTRACT 

This  report  describes  an  aircraft  trajectory  optimisation 
computer  program  that  has  been  used  successfully  at  Grunman 
for  a variety  of  problems  over  a period  of  many  years.  Three 
airplanes  are  simulated:  the  F-14,  F-15,  and  ATF  (Advanced 
Tactical  Fighter).  The  trajectories  are  in  three  dimensions 
over  a flat  earth.  Optimisation  is  achieved  by  means  of  the 
conjugate  gradient  variational  technique.  Inequality  constraints 
on  the  path  and  equality  constraints  on  the  final  point  are 
satisfied  by  penalty  integrals  and  penalty  functions,  respectively. 
The  equality  constraints  can  be  satisfied  to  the  limit  of 
computer  accuracy  by  a new  technique  that  avoids  Increases  in 
the  penalty  function  constants. 
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NOTATION 


constants  used  in  the  inequality  imposed  on  the  path 
by  the  structural  limit 

coefficients  varying  with  Mach  used  in  the  curve 
fit  for  the  lift  coefficient 

coefficients  varying  with  Mach  used  in  the  curve 
fit  for  the  drag  coefficient  used  when 
M>0. 9 

coefficients  varying  with  Mach  used  in  the  curve 
fit  for  the  drag  coefficient  used  when 
M s 0.9 

drag  coefficient;  see  Eqs.  (11)  and  (12) 
lift  coefficient;  see  Eq.  (10) 
function  of  Mach  appearing  in  Eq.  (12b) 
drag;  see  Eq.  (6)  and  Fig.  1 

time  derivative  of  the  state  vector  X;  see  Eq.  (21) 
acceleration  of  gravity 

equations  for  afterburner  blow-out  and  engine 
operational  limit  are  expressed  generically  by 

g (M,z)  - 0 

inequalities,  Eqs.  (7),  (8),  and  (9),  imposed  on  the 
path  and  the  mixed  state-control  inequality,  Eq.  (4), 
are  expressed  generically  by  % (X,u)  - 0 

Heavyside  unit  step  function  h (•)  is  zero  for 
negative  arguments  and  one  for  positive  arguments 

generalized  Hamiltonian;  see  Eq.  (27) 

the  true  Hamiltonian,  which  is  obtained  from  the 
generalized  Hamiltonian  by  setting  u to  the  value 
that  maximizes  the  latter  function 
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x,  y,  z 


diagonal  matrix  of  positive  constants;  see  Eq.  (36b) 
lift;  see  Eq.  (5)  and  Fig.  1 
aircraft  mass 

constant  appearing  in  the  limit  on  the  normal 
aerodynamic  force;  see  Eq.  (4) 

Mach;  see  Eq.  (7) 

limit  Imposed  on  Mach;  see  Eq.  (7) 

constant  appearing  in  the  limit  on  the  normal 
aerodynamic  force;  see  Eq.  (4) 

cost  functional;  see  Eq.  (23) 

cost  functional  for  a problem  with  constraints  on  the 
final  point 

augmented  with  penalty  terms;  see  Eq.  (36) 

cost  functional  for  the  fixed  final  time  approximation 
to  a problem  with  constraints  on  the  final  point; 
see  Eq.  (42) 

aircraft  reference  area 

time 

thrust 

control  vector  composed  of  a,  u,  and  q 
velocity  magnitude 

speed  of  sound;  see  Eqs.  (16),  (18),  and  (20) 

rectangular  coordinates  of  the  aircraft,  z is  the 
altitude 

vector  of  state  variables 

Initial  and  final  values  of  state  vector  as  stated 
in  the  problem  formulation 
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angle  of  attack  (deg.);  see  Fig.  1 
a + € (deg) 

upper  limit  on  a (deg);  see  Eq.  (3) 

engine  mass  flow  varying  with  M,  Z,  and  q 

vector  of  coefficients  used  in  Eq.  (32b);  see  Eq.  (33) 

coefficient  appearing  in  Eq.  (32b);  see  Eq.  (33) 

flight  path  angle  measured  in  the  vertical  plane 
between  the  velocity  and  the  horisontal;  see  Fig.  1 

vector  of  positive  constants  appearing  in  Eq.  (53) 

see  Eq.  (12b) 

sum  of  the  ram.  Interference,  and  spillage  drags  in 
the  F-15  engine 

fixed  angle  between  thrust  and  fuselage  reference 
line  (deg);  see  Fig.  1 

throttle  coefficient  appearing  in  Eqs.  (Id)  through  (lg) 
normalized  time;  see  Eq.  (2 A) 

Lagrange  multipliers  associated  with  the  state  vector  X 

adjoint  variables  associated  with  the  state  vector  X 

adjoint  variable  associated  with  the  free  parameter  x 
regarded  as  a state  variable 

bank  angle;  the  angle  between  the  lift  vector  L and  the 
vertical  plane  containing  the  velocity  vector 

atmospheric  mass  density;  see  Eqs.  (15,  (17),  and  (19) 

variable  coefficient  that  scales  the  increments  in  u; 
see  (Eq.  (32) 

final  time 
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final  times  for  problems  A,  B,  and  C;  see  Table  1 


X heading  angle  - the  angle  between  the  x axis  and  the 

projection  of  the  velocity  vector  on  the  horizontal 
plane 

T 

( ) transpose  of  a vector 

(*)  derivative  of  the  variable  with  respect  to  t 

II  I lR  generalized  norm  of  a vector;  see  Eq.  (36b) 

5(  ) variation  with  the  independent  variable  fixed 

A(  ) variation  with  the  Independent  variable  not  fixed 

( >€  ( ) belongs  to 

V(  ) for  all  ( ) 
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1.  INTRODUCTION 


This  report  des»_.ibes  a computer  program  that  calculates 
optimal  three-dimensional  aircraft  trajectories  over  a flat  earth. 
The  cost  functional  and  terminal  constraints  can  be  defined  at 
will.  At  present  there  are  three  FORTRAN  decks  simulating  the 
F-14,  F"15,  and  ATF  (Advanced  Tactical  Fighter).  A fourth  program 
simulates  the  evasion  of  a missile  by  an  F-14. 

The  optimisation  is  accomplished  by  the  conjugate  gradient 
variational  algorithm  (Ref.  1).  Constraints  on  the  final  point 
are  satisfied  by  a new  technique  that  avoids  increasing  the 
penalty  function  constants  (Ref.  2).  As  with  other  optimization 
techniques,  there  is  a possibility  of  converging  to  a local 
rather  than  the  absolute  minimum.  Whenever  a gradient- type 
algorithm  is  used  the  results  should  be  checked  by  computing 
the  problem  twice  using  different  starting  trajectories.  This 
not  only  tests  for  local  minimums  but  also  reveals  whether  or 
not  the  "noise"  inherent  in  numerical  methods  has  prevented 
convergence. 

Numerical  experiments  on  the  rate  of  convergence  for  two- 
point  boundary  value  problems  are  described  in  Ref.  2.  Some 
optimal  F-14  maneuvers  are  presented  in  Ref.  3. 

Descriptions  of  Grumman  programs  for  paths  over  a spherical 
earth  may  be  found  in  Refs.  4,  5,  and  6.  The  method  of  Refs.  4 
and  5 optimizes  trajectories  in  a vertical  plane  by  means  of 
gradients  (rather  than  conjugate  gradients)  and  satisfies  con- 
straints on  the  final  point  by  increasing  the  penalty  function 
constants.  The  method  of  Ref.  6 optimizes  three-dimensional 
trajectories  by  converting  the  variational  problem  to  one  with 
an  ordinary  minimum. 


2.  SIMULATION  OF  THE  F-14  TURBOFAN  AIRCRAFT 


DYNAMICAL  EQUALITY  AND  INEQUALITY  CONSTRAINTS 

The  three-dimensional  flight  of  the  F-14  turbofan  aircraft 
is  governed  by  the  following  differential  equations  (Ref.  7,  p.  48, 
p.  50).  Aircraft  and  thrust  sideslip  angles  are  neglected  and 
rigid  body  dynamics  are  not  included.  The  Earth  is  assumed  to  be 
flat.  Forces  and  angles  are  shown  in  Fig.  1. 
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Fig.  1 Forces  and  Angles  in  Vertical  Plane  when  Bank 
Angle  Equals  Zero 


x * v cos  7 cos  x 
y « v cos  7 sin  x 
z — v sin  7 
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(la) 

(lb) 

(lc) 


v-i  (“n  T cos  a'  - 1/2  p v2SCD)  — g sin  y (Id) 

V - C°8VU  (tj  T sin  a'  + 1/2  p v2SCL)  - ^ cos  y (le) 

* " W’wa^y  (T1  T sin  a'  + 1/2  p v2SCL)  (If) 

A - -P  (lg) 


Figure  1 shows  that  the  flight  path  angle  7 is  undefined 
when  the  magnitude  v of  the  velocity  vector  is  zero.  This 
is  associated  with  the  vanishing  of  the  denominator  of  Eq.  (le). 
Thus  the  initial  v must  not  be  specified  as  exactly  zero  for 
problems  that  begin  at  take-off. 

The  heading  angle  x is  also  undefined  whenever  the  projection 

of  the  velocity  vector  on  the  (x,y)-plane  has  zero  magnitude 

(i.e.,  v cos  y-  0).  By  definition  for  three-dimensional  problems 

y ranges  from  -90°  to  +90°,  and  it  is  extremely  unlikely  that 

the  limiting  values  would  ever  occur  in  practice.  However,  for 

motion  restricted  to  a vertical  plane  it  is  convenient  to  fix  x 

and  u at  zero  and  allow  y to  have  any  value.  Thus  cos  y can 

pass  through  zero.  Whether  computer  overflow  occurs  when 

calculating  Eq.  (If)  depends  upon  the  time  intervals  used  for 

the  numerical  integration.  The  computer  program  avoids  this 

-20 

overflow  by  resetting  cos  y to  ± 10  whenever  |cos  y|  is 

-20 

originally  less  than  10  . This  does  not  affect  the  trajectory 

since,  as  mentioned  above,  x is  not  obtained  from  Eq.  (If). 

The  control  variables  are  the  bank  angle  u,  thrust  coefficient 
q,  and  angle  of  attack  a (which  is  an  argument  of  C.  (M,a)  and 

I L 

a » a + e).  The  latter  two  control  variables  must  satisfy  the 
inequalities 

0 * r\  s 1 

a £ a 

max 
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(2) 

(3) 


The  limit  on  the  normal  aerodynamic  force  imposes  a mixed  state* 


control  Inequality. 

L cos  a + D sin  a £ mgn  if  m * in  (4a) 

L cos  a + D sin  a s mgn  if  m > m (4b) 

Here  the  lift  L and  drag  D are 

L - 1/2  p v2  SCL  (5) 

D - 1/2  p v2  SCD  (6) 

and  the  constant  n is  usually  chosen  to  be  between  6 and  8.  Mach 
number  is  defined  and  limited  by 

M - v/v,  ‘ ►W  O) 

The  structural  limit  imposes  a second  state  inequality,  which 
is  modeled  as 

2 

M < a + be  + cs  (8) 

Finally,  the  altitude  must,  of  course,  be  above  sea  level. 

z > 0 (9) 

ENGINE  CHARACTERISTICS 


The  thrust  T • T(M,a)  is  obtained  either  from  the  augmented 
thrust  table,  the  military  thrust  table,  or  set  to  aero,  depending 
upon  the  position  of  the  aircraft  relative  to  the  afterburner 
blow-out  altitude  and  the  engine  operational  limit.  The  data  for 
both  altitudes  are  curve-fitted  as  parabolic  functions  of  Mach. 

The  fuel  flow  0 es  P(M,s,t))  is  nominally  trivarlate,  but 
at  present  has  been  coded  only  for  tj  • 1 and  only  for  afterburner 
thrust.  This  fuel  flow  has  been  used  in  approximation  for 
trajectories  that  have  short  periods  either  with  n<l  or  with  mili- 
tary thrust.  The  data  for  the  engine  characteristics  were  obtained 
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from  Refs.  8 end  9.  The  interpolation  procedure  is  described 
in  Section  9. 

AERODYNAMICS 

The  lift  coefficient  CL(M,o)  is  curve-fitted  using 

CL  " *1^  + bl<M)  a (10> 

The  coefficients  e^  end  ere  obtained  from  tables  using  linear 
interpolation. 

When  M > 0.9,  the  drag  coefficient  CD(M,CL)  is  modeled  as 

CD  - a2(M)  + b2(M)CL  + c2(M)  Cl2  (11) 

At  lower  values  of  M the  relation  is 

CD  - a3(M)  + b3(M)  (ACl)2  +c3(M)  (ACl)4  + d3(M)  (ACL)6  (12a) 

ACL  * cl  ' CL  <12b> 

o 

Since  the  Macb-sweep  program  was  operative  during  the  tabulation 
of  the  original  data  (Ref.  10),  there  is  no  need  for  a wing-sweep 
variable. 
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3.  MODIFICATIONS  REQUIRED  FOR  THE  ATF  AIRCRAFT 

This  section  presents  the  changes  made  in  the  F~I4  program 
to  simulate  the  ATF  (Advanced  Tactical  Fighter).  The  fuel  flow 
0(M,  z,  n)  has  been  tabulated  for  both  afterburner  and  military 
thrust.  In  each  case  there  are  nine  values  of  the  throttle  setting 
t).  The  interpolation  is  performed  by  the  specially  coded 
subroutine  BVI3,  which  calls  on  the  subroutine  BVISP  described 
in  Section  9.  Interpolation  is  linear  in  the  z and  t]  directions 
and  by  spline  in  the  M direction. 

The  data  for  thrust  and  fuel  flow  were  provided  by  C. 

Giannetto  of  Engineering's  Navigation,  Guidance,  and  Control 

Section.  The  values  of  both  T and  p must  be  multiplied  by 

three  because  the  ATF  is  equipped  with  this  number  of  engines. 

These  data  must  also  be  multiplied  by  the  pressure  ratio 
2 2 

p(z)  v (z)/p(o)v  (0).  In  addition,  the  fuel  flow  data  are 

8 S 

multiplied  by  v (z)/v  (0).  The  engine  operational  limit  is 
s s 

60,000  ft  for  both  military  and  afterburner  power. 

The  formula  for  CL  differs  from  Eq.  (10)  only  in  that 
a^(M)  • 0,  i.e.,  CL  - b^(M)a.  The  equation  for  C^  has  the  form 
of  Eq.  (12)  with  C^  (M)  « d^(M)  a 0 and  with  the  store  drag 
ACp  added.  ° 

CD  " a3(M)  + b3<M>CL  + c3(M)CL  + ACd(M)  (13) 

This  equation  is  used  for  all  values  of  M. 

The  upper  limit  a on  the  angle  of  attack  is  no  longer 

max 

constant,  but  is  tabulated  as  a function  of  Mach.  As  presently 
coded,  the  trajectories  begin  with  military  thrust  and  light  the 
afterburner  at  a later  time.  The  latter  is  given  by  t - 
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r 


m 


where  d is  a predetermined  constant  and  t is  the  final  time. 

The  adjoint  variables  (see  Section  7)  do  not  jump  on  encountering 
a discontinuity  in  the  system  equations  that  are  defined  in  this 
way. 


4.  MODIFICATIONS  REQUIRED  FOR  THE  F-15 


Although  both  the  F-14  and  the  McDonnell- Doug las  F-15  use 
turbofan  engines,  the  ram,  interference,  and  spillage  drags 
cannot  be  neglected  when  simulating  the  latter  airplane.  The 
sum  of  these  drags  has  been  tabulated  as  AD(M,z).  Net  thrust 
is  obtained  by  a vectorial  subtraction  of  AD  from  gross  thrust. 
Since  AD  has  the  same  direction  as  D (Fig.  1),  Eq.  (Id)  is  replaced 
with 


$ - (qT  cos  a' -l/2pv2SCD  - AD) An  -g  sin  y (14) 

Equations  (le)  and  (If)  are  unchanged  in  form,  although  T now 
represents  gross  thrust. 

The  present  program  has  afterburner  thrust  data  only  for 
^ 1 and  has  no  data  at  all  for  military  thrust.  Thus  the  engine 

operational  limit  coincides  with  afterburner  blow-out.  This 
altitude  is  expressed  either  as  a linear  function  of  Mach  or  as 
75,000  ft,  whichever  is  lower.  The  equation  for  the  drag  coefficient 
Cp  has  the  form  of  Eq.  (12)  with  d^(M)  » 0.  The  same  equation  is 
used  for  all  values  of  M.  The  upper  limit  on  a is  tabulated  as 
a function  of  Mach, 
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5.  F- 14  (EVADER)  VERSUS  A MISSILE  (PURSUER) 


The  F-14  program  has  been  combined  with  a missile  program 
(Ref.  11)  to  simulate  missile  evasion.  The  control  variables 
of  the  F~14  are  chosen  so  that  it  attempts  to  evade  a pursuing 
missile.  The  latter's  strictly  deterministic  trajectory  is 
governed  by  proportional  navigation.  The  final  time  associated 
with  the  starting  trajectories  is  chosen  to  be  so  small  that  the 
F-14  cannot  be  overtaken.  At  this  final  time,  which  is  held 
fixed,  the  distance  separating  the  vehicles  is  maximized.  The 
minimum  separation  as  the  vehicles  move  along  the  final 
trajectories  is  found  using  parabolic  interpolation/ extrapolation. 
The  final  time  is  then  given  a small  increment  toward  the  time 
of  this  minimum,  and  the  process  is  repeated.  The  iteration 
converges  to  the  trajectories  that  either  maximize  the  minimum 
separation  or  maximize  the  time  of  impact. 

As  of  the  date  of  this  report,  the  full  simulation  of  the 
missile  trajectory  as  given  in  Ref.  11  has  not  been  completed  for 
the  above  program. 


6.  ATMOSPHERE 


Three  regimes  are  used  to  approximate  the  data  of  Ref.  12 

3 

for  the  mass  density  p (slugs/ft  ) and  the  speed  of  sound  v 

s 

(ft/sec).  The  following  equations  are  used  when  the  altitude 
is  below  36,146  ft. 

P-  2. 37688- 10_3(l-6. 7911- 10~6z)4, 3085  (15) 

v - 1116.45(1-6.863956- 10‘6t) 1/2  (16) 

s 

For  altitudes  between  36,146  and  65,874  ft  the  expressions 
are 

P - 3. 9792633- 10"3  exp (-4. 7829648- 10‘5z ) (17) 

vg-  968.08  (18) 

The  regime  for  higher  altitudes  is  modeled  as 

P-  5. 1526166*  10~3(1+1. 6606526- 10“6z)"32,838989  (19) 

vs-  922. 5793652(1+1. 535633914* 10“6z) 1/2  (20) 


7.  OPTIMIZATION  PROCEDURE 


The  choice  of  an  optimization  algorithm  was  largely  determined 
by  the  presence  of  the  tabulated  functions.  This  ruled  out 
second  order  methods  (e.g.,  second  variation,  quasilinearization) 
because  the  numerical  approximations  to  the  various  second  order 
partial  derivatives  are  unreliable.  The  classical  indirect 
and  the  MIN  H methods  could  not  be  used  because  it  is  not  possible 
to  express  explicitly  the  value  of  a that  maximizes  the  generalized 
Hamiltonian  (defined  in  Eq.  (27)).  To  avoid  these  difficulties 
the  gradient  method  was  tried  at  first.  Later  the  extension 
to  conjugate  gradient  (Ref.  1)  was  made  in  an  attempt  to  reduce 
computer  time.  A reduction  by  a factor  of  about  five  was  obtained. 

CONJUGATE  GRADIENT  ALGORITHM 

This  subsection  presents  the  conjugate  gradient  procedure 

for  problems  whose  final  point  and  final  time  can  be  varied 

freely.  On  defining  the  state  vector  h (x,y ,z,v, Y^m)  and 

T 

the  control  vector  u s (a,  p.,  q)  the  dynamic  equations 
(Eq.  (1))  can  be  written  as 

X - f(X,u)  (21) 

This  system  is  to  be  transferred  from  the  initial  point 

X(o)  - XD  (22) 

to  the  final  time  x and  final  point  X(x)  that  minimize  the  cost 
functional 

P - P(X(x),x)  (23) 

An  important  example  of  this  subsection' s problem  is  that  of 
maximizing  the  altitude  (P  - -z)  with  the  final  values  of  the 
other  components  of  X as  well  as  the  final  time  open. 


In  order  to  keep  u(t)  as  smooth  as  possible,  it  is  advisable 
to  avoid  extrapolation  by  replacing  t with 

e - t/x  (24) 

whose  range  is  always  0 to  1.  The  final  time  x could  now  be 
regarded  as  a free  parameter;  however,  the  analysis  will  be  a 
little  clearer  if  it  is  Instead  regarded  as  an  additional  state 
variable.  If  X continues  to  represent  only  the  original  state 
variables,  the  system  equation  takes  on  the  form 

dX/de  - x f (X,u)  (25a) 

dx/de  - 0 (25b) 

Note  that  x is  distinguished  from  the  other  state  variables  in 
that  its  Initial  value  is  not  subjected  to  a constraint  of  the 
Eq.  (22)  type. 

The  conjugate  gradient  algorithm  employs  adjoint  variables 
* , * that  obey 

A T 

~h  (Xx'>*,VT)  - x >£  f* 6u  ■ IS  6u  <26> 

The  second  member  has  been  simplified  by  defining  the  generalized 
Hamiltonian  as 

u 

H - T AX  f (X»U)  (27) 

The  differential  equations  that  * , A must  satisfy  when 

xx  J 

Eq.  (26)  holds  can  be  found  by  expanding  the  left  member  of  this 
equation  and  using  the  equations  of  variation  of  Eq.  (25). 

J 
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d*x/de  - - (dH/dX)T  (28a) 

d*T/de  - - 8h/St  (28b) 

When  the  final  values  of  the  adjoint  variables  are  defined  as 

*x(l)  - ' (dP/dX)T  , *t(1)-  - dP/dr  (29) 

the  Integral  of  Eq.  (26)  becomes 

-5P  - XT8X  + X Bt  + 

X T 

6mQ 

ri 

-xt(0)6t+  / |J  Oude  (30) 

The  last  member  assumes 

6X(0)  - 0 (31) 

since  all  the  trajectories  obey  the  initial  conditions  of  Eq.  (22). 
Equation  (30)  indicates  that  the  increments 


T 

6t  - oXT( 0)  , 5u(e)  - a( !g  (0))  (32a) 

are  in  the  direction  of  steepest  descent.  The  best  distance  to 
move  in  this  direction  can  be  found  from  a search  with  a. 

Although  the  first  descent  must  be  in  this  direction,  the  later 
descents  can  be  in  the  conjugate  gradient  direction  given  by 
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aj 


(0)  + eT 


(32b) 


Here  the  subscript  j represents  the  descent  number  and  and 
stand  for 


f>  - \2(o)/xJ  (0)  . a - 

j j-i 


/K-rn 


■1,  i 2 

l|B)  d8  03) 


Although  the  conjugate  gradient  theory  for  problems  with  n state 
variables  calls  for  recycling  after  n descents,  experiments  with 
the  present  variational  problem  indicate  that  Eqs.  (32a)  and 
(32b)  should  be  alternated  from  descent  to  descent  for  best  results. 


Whenever  Eq.  (32)  makes  a (0)  + 5a(0)  > a the  new  a 

max 

is  reset  to  a in  accordance  with  Ineq.  (3).  Similar  steps 
keep  tj  between  0 and  1. 

The  conjugate  gradient  algorithm  consists  of  the  following 
steps. 

(a)  Guess  u(©)  and  t 

(b)  Compute  X(©)  by  integrating  Eq.  (25).  Find  PQ 
from  Eq.  (23) 

(c)  Compute  *x(e)>  \(e)  by  integrating  Eq.  (28) 
backward  from  Eq.  (29) 

(d)  Obtain  four  values  of  a using  Eq.  (32a)  with  © «■  0 
while  5a,  6u , 6q,  6t  separately  take  on  small  values. 
Set  equal  to  the  smallest  of  the  a' s 

(e)  Increment  u(©)  and  t using  Eq.  (32)  with 
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(f)  Compute  X(6)  by  integrating  Eq.  (25).  Find 
P^  from  Eq.  (23) 

(g)  If  P^  > PQ,  a^/10  -*■  ct^.  If  the  new  o^  is  above  a 
very  small  tolerance  return  to  step  (e) , otherwise 
terminate  the  computation 

(h)  Set  02  m 4a^  and  find  P2  as  in  steps  (e)  and  (f) 

(i)  Set  a to  the  value  at  the  minimum  of  the  parabola 
through  (0,Po),  (av P^ , (a2,P2) 

f 0,1 f (42-l)P  -42p  +p  1 


' O f (4^-l)Po‘4ZP1-t-P2 

. 2 JL  (4“1)*0"4P1+  p2 


(j)  Obtain  a set  of  four  values  of  a as  in  step  (d) 
but  using  large  though  acceptable  values  of  6a, 

6u,  6t|,  and  5t.  Set  to  the  smallest  of  these 
and  the  a of  step  (i).  Compute  P^  as  in  steps  (e) 
and  (f) 

(k)  Set  PQ  to  the  smallest  of  P^,  P2,  P3«  Use  the 
corresponding  a in  Eq.  (32)  to  update  u,x.  If 
PQ  t P^,  recompute  X(0)  so  that  the  new  nominal 
trajectory  is  in  storage 

(l)  Return  to  step  (c) 

TWO- POINT  BOUNDARY  VALUE  PROBLEMS 

This  subsection  shows  that  the  solution  to  a problem  that 
specifies  constraints  on  the  final  point  can  be  approached  by  a 
sequence  of  variable  endpoint  problems  (Ref.  2).  First  the  time' 
optimal  case  will  be  treated. 


L 
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The  original  problem  --  to  be  called  Problem  A (see  Table 
1)  --  is  to  transfer  the  system  obeying  the  differential  equations, 
Eq.  (25),  and  the  initial  conditions,  Eq.  (22),  to 

X - Xf  (34) 

via  the  path  that  minimizes  the  cost  functional 

P*  - t (35) 

A 

The  final  time  defined  by  the  solution  extremal  will  be  called  r . 


TABLE  1 PROBLEM  FORMULATIONS 


Problem 

Cost 

Functional 

Initial 

Conditions 

Final 

Conditions 

A 

T 

X - XQ,  t open 
*x  open,  » 0 

X - Xf 

Xx  open,  *T  - -1 

B 

T + £|  IX  - X I I2 

K 

same  as 
above 

\ - -K(X  - X£) 

-1 

C 

lx  - xfl  I2 

K 

X - XQ,  x fixed 
\ open 

t open,  At»  0 
\ - -K(X  -Xf) 

The  algorithm  begins  with  a preliminary  stage  whose  purpose 
is  merely  to  find  an  optimal  trajectory  that  is  fairly  close  to  the 
solution  and  that  has  a lower  final  time.  This  is  accomplished 
by  solving  Problem  B (Table  1),  which  differs  from  Problem  A 
in  that  the  final  constraints,  Eq.  (34),  are  ignored  and  the 
cost  functional,  Eq.  (35),  is  replaced  with 

PB  - r + £1 IX(l)  - Xf I I2  (36a) 

K 
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The  penalty  term  Is  defined  as 


I IX-Xf I I2  = (X-Xf)TK(X-Xf)  (36b) 

X 

where  K is  a diagonal  matrix  of  positive  constants. 

After  Problem  B has  been  solved  by  the  procedure  of  the 

previous  subsection,  the  results  can  be  used  to  obtain  a first 

A 

order  estimate  of  t . Equation  (26)  is  integrated  assuming  only 

dH(0)/du  - 0. 


A 6X4 A 5t 

X T 


- X 6X4  A 6t 

X T 


e-o 


0-1 


(37) 


This  Important  relation  holds  between  the  endpoints  of  any 
extremal  and  those  of  any  adjacent  arc  whose  6X(0),  6x,  5u(0) 
obey  the  equations  of  variation  of  Eq.  (25).  Let  us  evaluate  it 
for  the  extremal  of  Problem  B using  the  endpoints  of  the 
(otherwise  unknown)  extremal  of  Problem  A.  This  Implies 

5X(0)  - 0 , 5X(1)  - Xf  - XB(1)  (38) 

so  that  Eq.  (37)  can  be  written  as 

6t  - X*(l)  [xf-3?(l)]  / [\(0)  - \<1)]  (39) 

In  accordance  with  Eqs.  (29)  and  (36) 

X*(l)  - - [X(l)  * Xf]TK  (40) 

so  that  the  first  order  estimate  for  the  difference  between  the 
final  times  of  the  two  extremals  is 
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(41) 


6t  - [xB(l)-XfJ  K [jCB(l)-XfJ  / [xB<0)-XB(l)j 

It  can  be  shown  with  the  aid  of  Eq.  (59)  that  5x  Is  positive. 

The  algorithm  now  proceeds  to  the  second  stage  which  works 
with  Problem  C (Table  1).  The  cost  functional 

PC  - 1/21 IX(1)  - Xfl I2  (42) 

K 

is  to  be  minimized  with  t kept  fixed  at  the  following  first-order 
A 

estimate  for  t . 


x^  » xB  + 6t  (43) 

The  solution  procedure  of  the  Conjugate  Gradient  Algorithm 
subsection  is  used,  except  that  6t  » 0 replaces  the  first  members 
of  Eqs.  (32a)  and  (32b)  and  the  second  member  of  Eq.  (29) 
requires  XT(1)  - 0 (see  Table  1). 

The  above  formulation  of  the  variable  endpoint  Problem  C 

has  been  chosen  so  that  its  solution  is  identical  to  that  of 

C AC 

Problem  A,  provided  t is  exactly  equal  to  x . If  t happens 

A 

to  be  smaller  than  t , the  boundary  value  errors  obtained  from 

the  solution  to  Problem  C may  exceed  the  tolerances.  A new 

A C C 

estimate  for  x can  then  be  obtained  from  x + 8x  -*■  x with  6x 

determined  by  Eq.  (41)  using  XC(1),  XC(0),  ^(1)  instead  of  XB(1), 

x\<0).  x®u>. 

r 

The  x defined  by  Eq.  (43)  might  also  be  too  large.  In 

this  case  the  algorithm  converges  to  a nonextremal  trajectory 

that  makes  the  boundary  value  errors  negligible.  There  is  then 

C A 

no  way  of  estimating  the  magnitude  of  x -x  . Therefore  the  following 
algorithm  recommends  that  the  increment  in  x given  by  Eq.  (41) 
be  divided  by  two  as  a safety  measure. 
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The  algorithm  is  summarized  in  the  following  steps: 

(a)  Compute  a starting  trajectory.  Set  the  descent 
counter  j to  zero 

(b)  Compute  a descent  cycle  for  Problem  B.  J + 1 -*  j 

(c)  Compute  Bx  from  Eq.  (41) 

(d)  If  j <5,  go  to  step  (b) 

(e)  If  the  conjugate  gradient  direction  decreases  x 
(see  Eq.  (32)),  go  to  step  (b) 

(f)  If  the  values  of  t + 6x  computed  from  the  three  most 
recent  descents  do  not  agree  to  within  a given 
tolerance,  go  to  step  (b) 

(g)  Compute  a new  starting  trajectory  with 
x + 6x/2  -*  r.  j «■  0 

(h)  Compute  a descent  cycle  for  Problem  C.  j + 1 j 

(i)  Terminate  the  computation  if  all  boundary  value  errors 
lie  within  given  tolerances 

(j)  Compute  8x-  from  Eq.  (41) 

(k)  If  j < 5,  go  to  step  (h) 

(l)  If  6x  < 0,  go  to  step  (h) 

(m)  If  the  values  of  Bx  computed  from  the  three  most 

recent  descents  do  not  agree  to  within  a given  tolerance, 
go  to  step  (h) 

(n)  Go  to  step  (g) 

The  numbers  five  and  three  in  steps  (d),  (f),  (k),  and  (m) 
were  chosen  arbitrarily.  Step  (e)  was  inserted  to  help  the 
program  converge  to  a xB  that  is  below  xA.  Step  (1)  helps  test 
for  convergence  because  Bx  approaches  a positive  value  when 


Literal  subscripts  indicate  partial  differentiation  in  this 
subsection.  The  independent  variable  is  once  more  t rather  than  6, 


C A C A 

t < t and  a negligible  value  when  t > t . 

These  methods  are  easily  extended  to  problems  other  than 
that  of  minimum  time.  For  example,  suppose  the  cost  functional 
la  X^l)  with  the  final  values  of  the  other  state  variables 
specified.  The  final  time  can  be  either  fixed  or  open.  Let 
us  define  X^  as  the  value  of  X^l)  on  the  solution  extremal.  If 
were  known,  the  gradient  program  could  generate  the  solution 
extremal  by  minimizing  the  error  norm,  Eq.  (32b)  with  X_  - X^. 

A n 

But  a first-order  approximation  for  XT  can  be  found  from  Eq.  (37) 
using  the  adjoint  variables  defined  by  the  extremal  that 
minimizes 


-X„(l)+l/2  l [VU-Xj  ]' 
i-1  L f 1 


(44) 


Numerical  results  for  F-14  trajectories  are  presented  in 
Refs.  2 and  3. 

DISCONTINUITIES  IN  THE  SYSTEM  EQUATIONS 


If  at  a time  t the  afterburner  blows  out  or  the  engine  reaches 
its  operational  limit,  X jumps  from  f”  - f(t”)  to  f+  - f(t+). 

If  the  nominal  arc  were  an  extremal,  the  vector  equation  that 
governs  the  jumps  in  the  Lagrange  multipliers  X would  be  well 
known  (Ref.  13,  p.  104).* 


X + 


«x  f+> 


(t$) 


Here  ft(X)  ■ 0 is  the  hypersurface  at  which  the  system  equations, 
Eq.  (21),  are  discontinuous  (see  Engine  Characteristics  in  Section 
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2).  Equation  (45)  Is  nonlinear  In  since  f+  depends  on  u+ 
which  for  extremals  is  a function  of  A+.  Taking  the  scalar 

X m • T 

product  of  both  sides  with  f shows  that  the  Hamiltonian  H «■  X f 
is  continuous. 

The  present  problem  is  slightly  different.  The  proper 
behavior  must  be  determined  for  adjoint  variables  A as  they 
are  integrated  along  arcs  with  Hu  ^ 0 and  u previously  defined. 

The  theory  of  the  gradient  procedure  requires  the  adjoint  variables 
to  jump  so  as  to  preserve  the  influence  function  properties 
expressed  by  the  following  equations  (cf.  Eq.  (30)). 

P(t)  - (dP(T)/dX(t'))  BX(t')  - - AT(t')BX(t')  , t'e  [0,t)  (46a) 

6u(t)  - 0,  Vt  € £•  ,t]  (46b) 

At  the  hypersurface  of  discontinuity 

AX  - 6X+  + f+6t  - 5X”  + f~6t  (47) 

Here  Bt  is  the  difference  in  the  times  at  which  the  nominal 
and  varied  arcs  reach  the  hypersurface  g • 0 and  the  6X*  are 
evaluated  with  time  fixed.  Thus  AX  lies  in  the  hypersurface 
tangent  space. 


$xAX  - gx6x“  + gxf'  5t  - 0 
Eliminating  Bt  between  Eqs.  (47)  and  (48)  yields 


6X_  + 


(48) 


(49) 


25 


The  definition  of  * must  indicate  that  any  OX  has  the  same 
effect  on  P as  the  corresponding  6X+. 


The  last  two  members  will  balance  for  an  arbitrary  OX  if 


Here  we  have  solved  for  \ rather  than  * because  adjoint 
variables  are  integrated  backward  in  time.  Note  that  the  right 
side  of  Eq.  (51)  is  Independent  of  ^ and  that  the  form  of  this 
equation  is  equivalent  to  that  of  Eq.  (45)  for  Lagrange  multipliers 

The  jump  in  A produces  a jump  in  the  Ou  defined  by  Eq.  (32). 
Thus  the  u of  the  converged  trajectory  should  have  a jump 
Identical  to  that  that  would  be  obtained  if  the  indirect  method 
were  used  instead  of  the  conjugate  gradient.  Note  that  the 
adjoint  variables  are  continuous  at  the  boundaries  of  the 
atmospheric  regimes  because  f is  then  continuous  (see  Section  6). 


INEQUALITIES  IMPOSED  ON  THE  PATH 


The  inequalities,  Ineqs.  (7),  (8),  and  (9),  imposed  on  the 
state  variables  will  be  satisfied  using  penalty  Integrals  (Refs. 
4 and  5).  This  method  is  readily  Incorporated  into  a gradient 
program;  hcwever,  it  provides  only  an  approximation  to  the 
optimal  path.  The  mixed  state-control  lnequalltly,  Ineq.  (4), 
will  also  be  treated  by  this  method  because  it  cannot  be  solved 
explicitly  for  the  control  a.  These  inequalities  can  be  written 


generic ally  as 


gj(X,u)  ^ 0 j • 1, ...  ,4  (52) 

The  system  equations  are  augmented  with  penalty  state  variables 

•ft 

that  obey 

dX.  9 

m 2 Mgj-ttj),  X±(o)  - 0,  j - i“8,  i-9 12  (53) 

Here  h is  the  Heaviside  unit  step  function,  and  the  positive 
constants  Aj  have  been  introduced  so  that  the  penalty  integrals 
begin  accumulating  tthile  the  trajectory  is  still  in  the  admissible 
region. 

A given  problem  is  solved  at  first  with  the  constraints,  Ineq. 
(52),  ignored;  that  is,  with  the  corresponding  penalty  function 
constants  in  the  cost  functional,  Eq.  (36),  set  to  zero.  If  the 
solution  violates  the  constraint  gj , a final  value  is  assigned  to 
x^(l)  (i  - j + 8).  This  value  is  of  course  smaller  than  that 
computed  for  the  unconstrained  trajectory.  The  modified  problem 
is  then  solved  with  the  new  boundary  condition  satisfied  by  the 
procedure  given  earlier  in  this  section  under  Two -Point  Boundary 
Value  Problems.  If  the  constraint,  Ineq.  (52),  is  still  violated, 
the  problem  must  be  solved  again  either  with  a smaller  final  value 
assigned  to  the  penalty  state  variable  X^l),  or  with  the  Aj_g 
made  more  positive,  or  with  a larger  penalty  function  constant. 

Note  that  the  right  sides  of  Eq.  (53)  are  continuous  even 
though  they  contain  step  functions,  so  that  by  the  theory  of 
the  preceeding  subsection  * is  also  continuous  when  • 0. 

The  analysis  (Ref.  14)  for  exactly  optimal  paths  also  requires  A 


The  original  components  of  the  state  vector  are  the  seven 
variables  of  Eq.  (1)  plus  t. 
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to  be  continuous  for  three  of  the  present  inequalities,  but  not 
at  the  junctions  between  interior  and  boundary  arcs  of  the 
"second  order"  constraint  z - 0.  But,  then,  the  penalty  integral 
technique  also  uses  a differing  differential  equation  for  A, 

Since  this  technique  employs  optimization  theory  to  bring  in 
the  constraints,  we  may  be  confident  that  it  yields  a good 
approximation  to  the  exact  minimum. 

INTERCEPTION  OF  A MOVING  TARGET 

This  subsection  presents  the  modifications  in  the  solution 
procedure  for  the  two “point  boundary  value  problem  that  are 
required  when  the  final  point  is  in  motion;  that  is,  when 
is  a given  function  of  t rather  than  a constant  (see  Problem 
A in  Table  1 earlier  in  this  section).  As  before,  the  original 
problem  will  be  solved  by  means  of  Problems  B and  C whose  cost 
functionals  now  have  the  arguments 

PB  = PB(X(l),Xf<T),T),  PC  * PC(X(l),Xf(x))  (54) 

For  Problem  C,  X(o)  and  t(o)  are  both  fixed  and  u(0)  is  to  be 

p 

chosen  so  that  X(l),  x(l)  minimize  P . Of  course,  we  know  that 
u(0)  has  no  effect  on  x(l),  but  transversality  conditions  are 
assigned  as  though  the  computer  does  not  know  this.  The  point 
is  of  no  significance  since  we  shall  see  that  the  relevant 
quantity  in  Problem  C is  Xt(o)  - XT(1),  and  this  is  independent 
of  XT(1).  Thus  for  both  Problems  B and  C 


The  second  of  Eq.  (38)  is  replaced  with 


so  that  Eq.  (39)  changes  to 


The  derivation  of  this  equation  does  make  use  of  the  relation 
5t(1)  • 8t(o). 


When  the  nominal  trajectory  has  converged  to  an  extremal 
with  a fixed  t (Problem  C) , it  can  be  shown  using  the  integral 
Eq.  (59).  of  the  next  subsection  that 


This  expression  is  valid  only  when  the  numerator  is  an  order  of 
magnitude  smaller  than  the  denominator.  However,  some  moving- 
target  problems  specify  some  of  the  components  of  f and  dX^/dr 
to  be  identical.  The  denominator  of  Eq.  (58)  then  has  a tendency 
to  be  small  so  that  Eq.  (57)  is  unreliable.  Such  problems  can 
be  treated  by  replacing  the  6t  of  Eq.  (57)  with  a fixed,  small 
positive  constant.  However,  this  means  that  the  optimal  final 
time  is  found  only  approximately  and  at  considerably  more  com- 
putational expense. 


FIRST  INTEGRALS 


The  state  variable  t is  of  course  constant  as  are  the  adjoint 
variables  associated  with  x,  y,  and  the  four  penalty  state 
variables  Introduced  above  in  the  subsection  Inequalities  Imposed 
cm  the  Path.  For  trajectories  that  have 
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converged  to  extremals,  the  Hamiltonian  H - t £ is  constant. 

i-1 

The  equation  d^/de  » -H/t  then  has  the  integral 


A + H6/x  » A_  + (A^  f )6  m const.  (59) 

Finally  the  quantity  A x - A y + A is  also  constant  along 

y x x 

extremals  (Ref.  15).  The  last  three  integrals  are  not  used  by 
the  computer  program  although  the  time  history  of  H is  printed 
in  order  to  verify  convergence. 


8.  NUMERICAL  INTEGRATION 
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The  differential  equations  are  integrated  numerically  using 
a modified  Adams  predictor-corrector  method.  The  procedure  is 
illustrated  for  the  vector  set  of  equations  X «■  f(X).  The  corrector 
equation  is 


*k+l  - Xk  + ¥ [* + f0tk+1>]  (60) 

There  are  two  corrections  at  each  time  interval.  For  the  first 
correction  Yfc+1  is  set  to  the  predicted  value 

Yk+1  mXk+f  &f<xk>  - £<*k-l>]  (61) 

At  the  initial  time  f(X_^)  is  assumed  to  be  equal  to  f(X  ).  For 
the  second  correction  Yk+1  is  set  to  the  Xfc+1  obtained  from  the 
first  correction.  Variations  in  the  time  interval  At  are  not 
permitted  by  the  present  program  because  values  from  the  forward 
integration  for  the  state  variables  are  required  during  the 
backward  integration  for  the  adjoint  variables. 
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9.  STORAGE  AND  INTERPOLATION  OF  TABULATED  FUNCTIONS 
AFTERBURNER  THRUST 

The  afterburner  thrust  T - T(M,z)  data  have  been  key-punched 
in  the  format  of  the  Grumman  Data  Systems  Software  Library 
subroutine  TAB IN  (Index:  12.7.0.2).  The  subroutine  itself 
is  not  used  to  read  in  the  cards  because  it  is  coded  in  single 
precision.  The  thrust  and  its  partial  derivatives  with  respect 
to  M and  z are  obtained  from  the  specially  coded  subroutine 
BVISP.  Spline  interpolation  with  respect  to  Mach  is  used  along 
the  curves  of  constant  altitude.  The  equations  are  in  the  form 
given  in  Ref.  16.  The  second  derivatives  at  the  mesh  points 
are  found  at  the  beginning  of  the  main  program  and  then  stored 
for  later  use  by  the  interpolation  subroutine.  If  the  thrust 
goes  to  zero  as  M decreases  to  some  value  called  M,  then  the 
second  derivatives  as  well  as  the  values  of  T are  set  equal  to 
zero  at  the  mesh  points  of  the  region  Me{o,M!]  . After  T has 
been  found  on  the  two  curves  that  straddle  the  aircraft' s 
altitude,  linear  interpolation  with  respect  to  z is  used  to 
obtain  the  thrust  of  the  vehicle. 

The  FORTRAN  statement  used  to  obtain  afterburner  thrust 
T(M, z)  is 

CALL  BVISP  (NO,  Nl,  N2,  Ll,  L2,  TAB,  M,z,  SP3,  DTDM,  DTDZ,  T) 

NO  is  the  fixed  point  number  that  has  been  assigned  to 
each  table.  For  afterburner  thurst  NO  is  1. 

Nl  is  the  block  of  storage  containing  the  number  of  first 
arguments  of  each  bivariate  table.  Nl(l)  contains 
the  number  of  Machs  tabulated  for  afterburner  thrust. 
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N2  is  the  block  of  storage  containing  the  number  of 
second  arguments  of  each  table.  N2(l)  contains 
the  number  of  altitudes  tabulated  for  afterburner  thrust. 

Ll  is  the  storage  block  that  contains  the  first  arguments 
of  each  table.  For  example,  Ll(l,2)  contains  the 
second  Mach  of  the  afterburner  thrust  table. 

L2  is  the  storage  block  that  contains  the  second 
arguments  of  each  table.  For  example,  L2(l,2) 
contains  the  second  altitude  of  the  afterburner 
thrust  table. 

TAB  is  the  storage  block  that  contains  the  values  of 
each  table.  For  example,  TAB (3, 2,1)  contains  the 
afterburner  thrust  corres ponding  to  the  third 
altitude  and  the  second  Mach. 

M is  the  current  value  of  the  first  argument. 

z is  the  current  value  of  the  second  argument. 

0 2 

SP3  is  the  storage  block  that  contains  3 T/3M  at  each 

mesh  point  of  the  afterburner  thrust  table.  For 

2 2 

example,  SP3(1,2)  contains  3 T/3M  at  the  second 
mesh  point  of  the  curve  for  the  first  altitude. 

DTDM  is  the  output  dT/dM. 

DTDZ  is  the  output  dT/dz. 

T is  the  output  afterburner  thrust. 

FUEL  FLOW  AND  MILITARY  THRUST 

The  fuel  flow  PCM,z,tj)  and  military  thrust  T(M,z)  have 


Recall  that  fuel  flow  has  been  key-punched  only  for  ry»l  except 
for  the  ATF.  Section  3 presents  the  interpolation  for  the  latter 
aircraft . 
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also  been  key-punched  in  the  format  of  subroutine  TABIN  (see  the 
previous  subsection.  Afterburner  Thrust).  Regions  for  which 
data  are  unavailable  are  indicated  by  zeros.  The  interpolation 
is  planar  and  is  computed  by  the  specially  coded  subroutine  BV1. 

If  values  are  defined  at  all  four  enclosing  mesh  points,  the  one 
that  is  furthest  from  the  nominal  point  is  discarded  and  the 
remaining  points  are  used  to  determine  the  interpolating  plane. 

In  the  other  cases,  the  three  closest  mesh  points  that  have  defined 
values  and  are  not  collinear  are  used. 

No  significant  changes  in  the  trajectories  were  observed 
when  afterburner  thrust  was  obtained  by  spline  rather  than 
the  original  planar  interpolation.  Since  the  trajectories  are 
even  less  sensitive  to  the  tables  of  this  subsection,  their 
interpolation  was  kept  planar. 

The  FORTRAN  statement  used  to  obtain  either  fuel  flow  or 
military  thrust  is 

CALL  BVI  (NO,  Nl,  N2 , Ll,  L2 , TAB,  M,  z,  DADM,  DADZ,  A) 

The  arguments  are  the  same  as  those  of  BVISP  except  that  SP3 
is  omitted.  NO  is  4 for  fuel  flow  and  5 for  military  thrust. 
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APPENDIX 


COMPUTER  PROGRAM  LISTING  AND  TYPICAL  SOLUTION 

The  FORTRAN  IV  Computer  program  listing  is  reproduced 
on  the  following  pages  with  all  classified  statements  and  data 
replaced  by  asterisks.  The  particular  problem  is  the  minimum 
time  to  climb  of  the  F-1A  from  the  end  of  the  runway  to  z (altitude) 
“ 55,000  feet  and  V (velocity)  ■ 15A8.8  feet/second,  with  z con- 
strained to  be  always  positive.  The  numerical  solution  to  this 
problem  follows  the  program  listing. 

For  those  interested  in  seeing  the  classified  version  of  the 
above  problem,  the  program  listing,  including  propulsion  and 
aerodynamic  tables,  is  available  from  the  author. 
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